home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byte1286.arc / AITKEN.BAS next >
Encoding:
BASIC Source File  |  1986-10-29  |  2.8 KB  |  87 lines
  1. 100 DEFDBL A-H,L,O-Z
  2. 110 REM Use EPS=1E-6 if functions are single precision.
  3. 120 EPS=1E-10
  4. 130 SQT3=SQR(3)
  5. 140 DIM TABLE(20)
  6. 150 LG10=LOG(10)
  7. 160 DEF FNL10(X)=LOG(X)/LG10
  8. 170 F$="f(x)=SQR(EXP(x)-1)/SIN(x)"
  9. 180 DEF FNF(X)=SQR(EXP(X)-1)/SIN(X)
  10. 190 T1$="######   #####.##########    ####    ---"
  11. 200 T2$="######   #####.##########    ####    ###"
  12. 210 T3$="######  #####.##########    ---"
  13. 220 T4$="######  #####.##########    ###"
  14. 230 PRINT "Approximating the integral of the function:"
  15. 240 PRINT "       "; F$
  16. 250 INPUT "Number of lines in Gauss column (>1)";KL
  17. 260 INPUT "Lower limit for integral";A
  18. 270 INPUT "Upper limit for integral";B
  19. 280 IF B<=A THEN PRINT "Error. Lower limit <= upper limit.": STOP
  20. 290 PRINT
  21. 300 GOSUB 410
  22. 310 PRINT
  23. 320 FOR J=1 TO 20
  24. 330 INPUT "Enter 1 for next Aitken column, 0 to quit";CQ
  25. 340 IF CQ=0 THEN J=20: GOTO 390
  26. 350 GOSUB 760
  27. 360 PRINT
  28. 370 KL=KL-2
  29. 380 IF KL<3 THEN J=20
  30. 390 NEXT J
  31. 400 END
  32. 410 IF KL<1 THEN PRINT "Error in Gauss subroutine: KL<1.": STOP
  33. 420 NLINES=KL
  34. 430 NSUBS=1
  35. 440 NFUNCT=0
  36. 450 PRINT "Gaussian integration table for the function"
  37. 460 PRINT F$; ",  for x= ";A; " to ";B;"."
  38. 470 PRINT "Table contains ";NLINES; " lines."
  39. 480 PRINT
  40. 490 FOR JLINE=1 TO NLINES
  41. 500 XH=(B-A)/NSUBS
  42. 510 XH2=XH/2
  43. 520 XR=XH2/SQT3
  44. 530 START1=A-XH2-XR
  45. 540 START2=A-XH2+XR
  46. 550 SUM=0
  47. 560 FOR K=1 TO NSUBS
  48. 570 SUM=SUM+FNF(START1+K*XH)+FNF(START2+K*XH)
  49. 580 NEXT K
  50. 590 SUM=SUM*XH2
  51. 600 NFUNCT=NFUNCT+2*NSUBS
  52. 610 IF JLINE>1 THEN 650
  53. 620 PRINT "     N       Gauss             NF  Est.S.D."
  54. 630 PRINT USING T1$; NSUBS,SUM,NFUNCT
  55. 640 GOTO 720
  56. 650 RELERR=EPS
  57. 660 IF SUM<>0 THEN RELERR=ABS(TABLE(JLINE-1)-SUM)/ABS(SUM)
  58. 670 IF SUM=0 AND TABLE(JLINE-1)<>0 THEN RELERR=ABS(TABLE(JLINE-1)-SUM)/ABS(TABLE(JLINE-1))
  59. 680 IF RELERR<=0 THEN RELERR=EPS
  60. 690 NUMSD=-FNL10(RELERR)
  61. 700 IF NUMSD<0 THEN NUMSD=0
  62. 710 PRINT USING T2$;NSUBS,SUM,NFUNCT,NUMSD
  63. 720 NSUBS=2*NSUBS
  64. 730 TABLE(JLINE)=SUM
  65. 740 NEXT JLINE
  66. 750 RETURN
  67. 760 IF KL<3 THEN PRINT "Error in Aitken: KL<3.": STOP
  68. 770 KLM2=KL-2
  69. 780 PRINT: PRINT "Aitken column. ";KLM2;" entries."
  70. 790 PRINT: PRINT "   Line     Aitken          Est.S.D.
  71. 800 KLM1=KL-1
  72. 810 FOR JLINE=2 TO KLM1
  73. 820 TOP=(TABLE(JLINE+1)-TABLE(JLINE))^2
  74. 830 BOT=TABLE(JLINE+1)-2*TABLE(JLINE)+TABLE(JLINE-1)
  75. 840 IF BOT=0 THEN TABLE(JLINE-1)=TABLE(JLINE+1) ELSE TABLE(JLINE-1)=TABLE(JLINE+1)-TOP/BOT
  76. 850 JLM1=JLINE-1
  77. 860 IF JLINE<=2 THEN PRINT USING T3$;JLM1,TABLE(JLM1): GOTO 940
  78. 870 RELERR=EPS
  79. 880 IF TABLE(JLM1)<>0 THEN RELERR=ABS(TABLE(JLM1-1)-TABLE(JLM1))/ABS(TABLE(JLM1))
  80. 890 IF TABLE(JLM1)=0 AND TABLE(JLM1-1)<>0 THEN RELERR=ABS(TABLE(JLM1-1)-TABLE(JLM1))/ABS(TABLE(JLM1-1))
  81. 900 IF RELERR<=0 THEN RELERR=EPS
  82. 910 NUMSD=-FNL10(RELERR)
  83. 920 IF NUMSD<0 THEN NUMSD=0
  84. 930 PRINT USING T4$;JLM1;TABLE(JLM1);NUMSD
  85. 940 NEXT JLINE
  86. 950 RETURN
  87.